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ABSTRACT 

Modified gravity models require a screening mechanism to be able to evade the stringent constraints from local gravity experiments 
and, at the same time, give rise to observable astrophysical and cosmological signatures. Such screened modified gravity models 
necessarily have dynamics determined by complex nonlinear equations that usually need to be solved on a model-by-model basis to 
produce predictions. This makes testing them a cumbersome process. In this paper, we investigate whether there is a common signature 
for all the different models that is suitable to testing them on cluster scales. To do this we propose an observable related to the fifth 
force, which can be observationally related to the ratio of dynamical-to-lensing mass of a halo, and then show that the predictions 
for this observable can be rescaled to a near universal form for a large class of modified gravity models. We demonstrate this using 
the Hu-Sawicki f(R), the Symmetron, the nDGP, and the Dilaton models, as well as unifying parametrizations. The universal form is 
determined by only three quantities: a strength, a mass, and a width parameter. We also show how these parameters can be derived 
from a specific theory. This self-similarity in the predictions can hopefully be used to search for signatures of modified gravity on 
cluster scales in a model-independent way. 

Key words, cosmology: large-scale structure of Universe - cosmology: dark energy - gravitation - galaxies: clusters: general - 
galaxies: kinematics and dynamics 


1. Introduction 


A possible solution for the observed acceleration of our Universe 
is that general relativity (GR) is not the correct theory to describe 
gravity on large scales. Over the past decade, a large number of 
alternative gravity theories have been proposed ( |Clifton et al.| 
2012) 1. Many of these models cannot alleviate the existing, or 
they introduce new problems; nevertheless, the idea that cosmic 
acceleration has a gravitational origin persists. 


High precision tests of gravity indicate that general relativ¬ 
ity is the correct theory for describing gravitational physics on 
Earth and in the solar system ( |Bertotti et al. 2003 1 Will 2006| 
Willi ams et al,|2004 1 . This places strong constraints on any the¬ 
ory that seeks to modify general relativity. Thus, if gravity is 
modified on large scales - giving rise to cosmic acceleration - 
then a way of suppressing the modifications of gravity in the 
well-tested regimes is required ( |Khoury||20131 >. A key feature 
of screening mechanisms currently being considered is that they 
are fundamentally nonlinear, thereby making the study of their 
cosmological e ffects challenging (|Vainshtein |1972| |Khoury &| 
Weltman|2004[|Hmterbichler & Khoury|2010 1 . 


As previously argued in Gronke et al. ( |2014| |2015| >, devia¬ 
tions from GR occurring in screened-modified gravity theories 
(SMGs) can generally be grouped into three categories: no devi¬ 
ations (fully screened regime), maximum deviations (unscreened 
regime ), and the intermediate deviations (partially screened 
regime). As a natural continuation of this idea, the main ques¬ 
tions we would like to address in this paper follow: 


- How do the extents and occurrences of the three regimes vary 
with the screened modified gravity model and its parameters? 

- Is it possible to remap observables within one model and/or 
between different models; i.e., are SMGs self-similar in some 
way? 


Previous work has been carried out in this direction. IBrax et a~Ll 
( 2012a|b l have developed a theoretical framework in which 
it is possible to describe certain SMGs with two free func¬ 
tions. Using this framework |Winther & Ferreira|f20T5| propose 
an approximate scheme to perform fast /V-body simulations of 
SMGs by combining linear theory with a screening factor ob¬ 
tained from studying spherical symmetric systems. Also, [Mead] 


et al. (2015) have recently published a technique for remapping 
ACDM power spectra to the Chameleon f(R) model with ~ 3% 
accuracy. Complementary to that, we want to base this work on 
an observationally focused approach. Specifically, in this paper 
we focus on the behavior of screened modified gravity theories 
inside and in the vicinity of clusters of galaxies. In general we 
expect the strongest signal of modified gravity to be found in 
cosmological voids, but clusters are easier to work with theo¬ 
retically and observationally, and they remain an important ob¬ 


servational probe of gravity on intermediate scales (Baker et al. 
|20l5]>. 


To answer the questions above, we have structured this paper 
as follows. In Sec.[2j we give a general introduction to the topog¬ 
raphy of SMGs. We also picked three models as examples, which 
we describe in more detail. This includes the respective solving 
method for halo density profiles. In Sec. [3] we present the nu¬ 
merical results of these three models. In addition, we study the 
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effects of varying the model parameters on the observables and 
try to empirically find a rescaling that counterbalances this vari¬ 
ation. Also in this section, we develop a general analytic rescal¬ 
ing method and compare it to our numerical results. In Sec. [4} 
we discuss our results and possible observational implications. 
We conclude in Sec.0 

Throughout this work we use (Q m o, Hq) = (0.3, 0.7), M p 2 = 
87tG, p c = 3H 2 M 2 V p/) - Q. m p c and denote values today with a 
subscript zero. 


2. Methods 

In this section we give a brief introduction to screened modified 
gravity theories. We present a few example models in more detail 
and then give a general framework for considering a large class 
of models. In the end we discuss the numerical implementation 
for solving for the field profiles of the models for the case of a 
Navarro-Frenk-White (NFW) density profile. 


2.1. Modified gravity theories 

When we speak about modification of gravity or modified grav¬ 
ity theories, we generally mean the addition of terms to the 
Einstein-Hilbert action. Specifically, in scalar-tensor theories 
one (or multiple) additional scalar field is added with their re¬ 
spective Lagrangian L v , and the general action reads as 


■/ 


S = d 4 x 


M 2 , ' 

—pR + L v {ip, dip, d~ip) 


+ S ffl 0A (!) ,^ v ), (1) 


where the matter fields i/j (n couple to the metric g MV = A 2 g MV , and 
R is the Ricci scalar. We limit ourselves to one scalar field and a 
conformal coupling to matteiQ Since the scalar field is coupled 
to matter, i.e., A = A(ip), an extra fifth force arises that is in the 
non-relativistic limit, and per unit mass is given by 


= -V log A . 


( 2 ) 


The predictions of general relativity have been confirmed to 
high accuracy in the laboratory and in the solar system, which 
severely constrains any new such force in our (Earth’s) local en¬ 
vironment ( |Bertotti etal.|2003)|Will|2006[ [Williams et al.|2 004). 
If the fifth force is weak close to earth, then one would naively 
expect it to also be weak on cosmological scales. However, in 
the past decade several theories have been proposed where one 
is able to dynamically suppress the effect of such a fifth force in 
high density environments (relative to the cosmic mean) and still 
have interesting cosmological signatures. If this is the case, then 
a theory is said to possess a screening mechanism (for reviews 
see, e.g., Clifton et al.|2012 Khoury |20 10 Joyce et al.|2014 1 . 

The action Eq. (jTJ allows a loose mathematical classification 
of possible screening mechanisms: 

1. Screening due to the scalar field value ip. If L^ possesses 
a potential V(ip), it is possible to suppress the fifth force 
thanks to the value of the scalar field moving inside the po¬ 
tential. Physically, this can occur if the coupling strength de¬ 
pends on the local environment as is the case in the Sym- 
metron mechanism ( Hinterbichler & Khoury j|2010| Hinter- 
|bichler et al. [2011) , Another possibility is that the range of 
the fifth force depends on the environment, which is often 


1 Disformal coupling allows for additional screening mechanisms 
i Koivisto et al.|20 1 2|. 
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Table 1 . Overview of screened modified gravity models 


Screening type A- body code ' 


Model 


f(R) HS (a) 

Dilaton 

Symmetron 

DGP 

Galileon 

k-mouflage 


(1) Chameleon 
(1) Chameleon-like 

(1) Chameleon-like 

(2) Vainshtein 
(2) Vainshtein 

(2) Vainshtein-like 


(E), (I), (MGG), (O) 
(E) 

(E), (I) 

(E), (S), (K) 

(E), (I) 


(a) |Hu & Sawicki|(|2007t f(R)~ 

(A) ( E) ECOSMOGlLi et al |2012| >, (I) ISIS ( |Llinares et al.|2014] 
(K) |Khoury & Wyman| (|2009|), (MGG) MG-GADGET (|Puchwein| 
et al.||2013D, (O) |Oyaizu| (|2008), (S) |Schmidt et al.| (|2009|l; 
Schmidt| ( |2QQ9] > 


called the Chameleon mechanism (Khoury & Weltman 2004] 


Khoury & Weltman 2004 |Mo ta & Shaw 2007], A frequently 
considered SMG possessing the Chameleon screening is the 
[Hu & Sawicki|(20O7| f(R) formulation. 

2. Screening due to derivatives of the scalar field dip and/or d 2 ip. 
In this case, the screening is caused by the relation between 
the kinetic term of L v , i.e., -j(dip) 2 , and the rest of the scalar 
field Lagrangian. A popular mechanisms incorporating this 
idea is the Vainshtein mechanism ( | Vainshtem]| 1972[ ), where 
screening depends on the value of d L ip, and the k-mouflage 
mechanism (Babichev et al. |2009]|Brax & Valageas|2014a| , 
where the magnitude of (dipY decides whether a region is 
screened or not. Particular SMGs employing the Vainshtein 
Mechanism are massive gravity ( de Rham|2014|l, Gableon s 
( |Nicolis et al.|2009] >, and the DGP model ( |Dvah et al.|2 000). 

A selective overview over some models and their employed 
screening mechanism is given in Table [I] 

Over the past decade there have been several studies of mod¬ 
ified gravity on nonlinear scales, often using A'-body simula¬ 
tions. For theories that screen via the first method described 


above (chameleon-like screening), there have been studies of 

f(R) gravity (Li et al.|2012 Llinares et al. 

2014 Puchwein et al. 

2013] Zhao et al. 2011a; |Oyaizu 2008 

le et al. 20141 Lom- 


the Symmetron ( Davis et al.||201 2] |Llinares et al. 2014] Brax| 


|et al.||2012c|>, and the environment-dependent dilation model 


(Brax et al. 20TT| 2012c|l. For models that employ the second 


screening method (Vainshtein-like screening), there have been 


studies of the DGP model ( 

Schmidt et al. 2009; 

Schmidt 2009; 

Chan & Scoccimarro 2009 

Khoury & Wyman 

2009 Li et al. 

2013b 

Falck et al. 2015 2014), Gableons (Barreira et al. 2013) 

Li et al.|20l3a 

), and k-moflage (Brax & Valageas 2014b 

). This 


is by no means a complete list. 

In spite of the differences between the different screening 
mechanisms and the theories employing them, some common 
quantities can be defined. In particular, we focus on the enhance¬ 
ment of the gravitational force 


y = Geff/G - 1 


(3) 


with its theoretical maximum value y max , given a specific modi¬ 
fied gravity model. As described in j ]4.2[ y can be related to the 
ratio of the dynamical to the lensing mass. 


2.2. Example models 

The screened modified gravity models that we considered pos¬ 
sess three intrinsically different screening mechanisms: the Vain¬ 
shtein (Vainshtein]| 1972] screening where the derivatives of the 
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scalar field play a major role; the Symmetron (Hinterbichler et al. 
2011 1 Hinterbichler & Khoury|20T0]l screening, where the cou¬ 


pling strength is environment dependent and goes to zero in high 
density regions; and the Chameleon (|Khoury & Weltman|2004[ 


Khoury & Weltman 2004) screening, where the range of the fifth 


force approaches zero for regions with high matter density. As 
an example of the latter, we focus on the |Hu & Sawicki| (2007) 
f(R) model, which incorporates the Chameleon screening mech¬ 
anism. 

In the following, we do not present the example models in 
great detail. Instead, we refer the reader to our previous work 
( Gronke et al.||2014 20151, where we stated the full equations 
for the Symmetron and the |Hu & Sawicki| ( |2007[ ) f(R) model 
or - even better - to excellent reviews, such as those by Clifton 


et al. 
( 2015 


( 2012 ) , 


|Khoury (2013l, Joyce et al. ( 2014| >, and, Koyama 


2.2.2. Hu-Sawicki f(R) model 

The f(R) gravity models are a group of modified gravity theories 
where the Ricci scalar R in the Einstein-Hilbert action is replaced 
by a generic function R + f(R). By applying a conformal trans¬ 
formation, the f(R) action can be brought in the form of Eq. Q 
(see, ej» |Clifton et al.||2012| |Brax et aI7||20081 >. This transfor¬ 
mation 2 ] translates /(R)-gravity into a scalar tensor theory. The 
theory screens via the chameleon mechanism. 

The fifth force is given by 

/V = -\vfn (12) 

where fa = d f(R)/ d R is the scalar field. In a FLRW background 
spacetime and in the quasi-static limit (Noller et al. 2014) , the 
field-equation becomes 


2.2.1. Symmetron model 

In the Symmetron model (Hinterbichler & Khoury |20T0 ; Hinter- 
bichler et al.|201 1 ) , both the coupling function and the potential 
are symmetric around cp - 0 are given by 


A( ") = 1 + ^ 
V(<p) = -\r 2 <p 2 + 


(4) 

(5) 


Here, the prefactors M, p and A can be rewritten as 


M 2 _ ^'^‘mOPcoL 


ssb 


P = 


V2L 


A 1/2 = 


a ssb M Pl 


V8L 3 G,„ 0 p c o^, 


( 6 ) 

(7) 

( 8 ) 


which leaves the Symmetron parameters to be the scale factor at 
symmetry breaking c/ ssh , the range of the fifth force in vacuum L, 
and the fiducial coupling j3. 

In a Friedmann-Lemaitre-Robertson-Walter (FLRW) back¬ 
ground and in the quasi-static limit (Noller et al.|2014[ >, the field 
equation becomes 


V 2 <p = 

2A 2 


~ssb Pm - -3 

= -f + f 

a 3 Pm 


(9) 


where ip = yAip/iJ.. The fifth force is given by 


/> = 


ipVip 

ML 


( 10 ) 


In regions of high ambient matter density (p,„ » p 2 M 2 = 
pco/a\), the field will reside close to cf> — 0, and the fifth force 
will be screened. In low-density regions, on the other hand, the 
fifth force can be in full operation. In the Symmetron model, the 
maximum enhancement of the gravitational force, y max , is given 
by 


7max = 2/J 2 




(ID 


Since /? is a free parameter, the model can, in principle, give rise 
to an arbitrarily large (or small) deviation. How the screening 
mechanism works is explained in more detail in Sec.|2.3| 


r2 2 / Pi 


1 +4- 


a, 


fm \ 

f* 


Pn 
n +1 


V z /« = -& m H- 0 a z I = - 11 + a 2 H 2 n m x 


■a- 3 -4 iiA 


I Ga 

a 


(13) 


where n and fan are dimensionless model parameters. The max¬ 
imum enhancement of the gravitational force, i.e. the enhance¬ 
ment when there is no screening, is simply 


= 2/J 2 = -, 


(14) 


so gravity can be enhanced by up to 33%. How the screening 


mechanism works is explained in more detail in Sec. 2.3 


2.2.3. DGP 

The |Dvali, Gabadadze, & Porratij ( |200Q} > (DGP) model is an ex¬ 
ample of a braneworld model where we are confined to live 
in a four-dimensional brane that itself is embedded in a five¬ 
dimensional spacetime. The DGP model can give rise to self¬ 
acceleration of the Universe without an explicit cosmological 
constant, but this branch has problems with instabilities. We take 
the normal branch (n)DGP model as our main example. To get 
accelerated expansion in this branch, we need to add dark energy, 
and for simplicity, we assume that the background expansion is 
the same as in ACDM. This choice is not expected to change any 
of our conclusions. 

Since it is a higher dimensional theory, the DGP model only 
fits into the formalism of Eq. (jTji as an effective theory. However, 
the modifications of gravity in this model can also be described 
as a fifth force given by 

F v = (15) 


where the evolution of the scalar field <p, the co-called brane- 
bending mode, is determined by the field equatior0 


vV 


3/?dgp (a)fl‘ 


((V^) 2 - (V,V # ) 2 ) = 


O H- 
a /?DGp( a ) 


(16) 


2 g)ir ~ 4 2 (i P)g li v where A(tp) = e^ /Mpl with /? = 1/ V6 

3 For the self-accelerating branch, the same expression holds but with 

r c —> —r c . In this case, H(a)/H 0 = + yj(il m a~ 3 + (4 with 
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where 


H(a) 


PoGp( a ) - 1 + 2 (r c H 0 )—— 1 + — — -r 


H 


H 0 \~ 3 H 2 

The maximum value of the gravitational force in this model is 
1 


(17) 


Tmax 


3/?DGP 


(18) 


For high values of r c , we have /?DGp( a = 1) oc r c , so y max oc 
The larger r c , the weaker the effect of the fifth force. 

As also stated in [Schmidt] (20101, in spherical symmetry we 


can integrate Eq. (]T6]> directly to obtain 




2c 


r 3yS D GP (a)a 2 \ r ) a/3 DCP (a) 


& m Hl Jo<5mr 2 dr 


(19) 


This translates to a solution for the fifth force profile given by 


r = 


V 


l 


2(— 1 + VTT7(7)) 


e(r) 


*3?' 3/i|)fip 

where 

8(r c H 0 ) 2 Cl m £ S m r 2 dr 


e(r) = 


3/3l Cp (a) a3 


( 20 ) 


( 21 ) 


We have two regimes. For e <s 1 (large r), we get y = 
and gravity is maximally modified. On the other hand, for e » 1 
(small r), we have 


1 


1 


r = 




3/?dgp \[e 3/5 D gp ’ 

and the fifth-force is screened. 


( 22 ) 


2.3. The unified { m(a)./3(a )) description of chameleon-like 
models 


Brax et al. ( |2012b|a) ) show that any scalar-tensor theory that 
screens using the Chameleon and/or Symmetron mechanism can 
be described uniquely by the evolution of the mass and matter 
coupling along the cosmological attractor. Thus parametrizing 
these two functions, which have a clear intuitive meaning - the 
range and the strength of the fifth force, respectively - to be ef¬ 
fectively parametrized. The mapping between the potential V(cp) 
and the coupling function A(ip) to {m(a),/3(a)} has been derived 
in Brax et al. (2012b| and reads as (in parametric form) 


V(a) =Vo - —y 

M-i 


P 2 (a) 


<p(a) =(fii + 


f 

J a , am 2 (a) 

r m 

J a; 


Pm(a) d a 


M P \ J a am 2 (a ) 


Pm(a ) da 


log A(a) = log A(a,) + — r 


f 

J di 


P~(a) 


M 2 , Ja, am 2 (a) 


p m (a)da, 


(23) 

(24) 

(25) 


where p m (a) = 3H 2 M 2 l Tl m /a i . 

One of the advantages of using this form is the direct re¬ 
lationship between the evolution of linear matter perturbations 
and the fi, m functions 


dm 2 Hd m — ^ Tl nl u 


dm 1 + 


2 p\d)k 2 
k 2 + a 2 m 2 {a) 


(26) 
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The two functions are, in the case of the Symmetron model, 
given by 


m(a) 

m 

i.e., they follow the same evolution. The parameters Pq and mo 
are related to the model parameters described in the previous 
section via Po - P an d mo = j-. 

For the Hu-Sawicki f(R) model, on the other hand, we have 

m - 4= (28 > 

V6 

and 



H 0 VQ,„ + 4Q a / Tl m a~ 3 + 4Q A \" /2+2 
o/Um\(n + 1) \ n,„+4fi A / 


(29) 


which means that p is constant, and m(a ) oc a - 3 ( ,! /2+i) f or sma ll a. 
Note that the nDGP model does not map on to this parametriza- 
tion except for the case of linear perturbations where it would 
correspond to m(a) = 0 and p(a ) = 1 / -\/6/5dgp- 

When the scalar field, in a region of density p is close to the 
minimum of the effective potential (this corresponds to the fully 
screened regime), we have tp ~ ip(a p ) where a p = a emr (pT,/p) 1 ^. 
The screening condition for an object with Newtonianpotential 
<I* v in a region of density p can therefore be written a>|j 


e(a p , fl e nv) = min 


I ‘P(cip) ~ p(a e nv)[ 1 

3-P(ucnv ) 4/p] 


(30) 


An object is screened whenever e(a p ,a env ) •« 1. Here, a env = 
«(P^/Penv) 1/3 where p env is the density of the environment the 
object is located in. The min condition ensures that we get the 
correct value, e — 1, in the nonscreened regime. The screening 
condition is related to y via 

y — 3.p~(u env )e(ci p , ci e n v ). ( 31 ) 


We should note that this is a very rough analytical approxima¬ 
tion, but it is usually good enough for order-of-magnitude es¬ 
timates. To get accurate predictions, we need to solve the field 
equation numerically, as discussed in the next section. 


2.4. Numerical solving methods 

To obtain accurate predictions for the field and force profiles in 
the models we study we need to numerically solve the field equa¬ 
tions Eqs. <[9]» and ( fl3| ). We solve these equations by discretizing 
them on a grid, using the NFW density profile for p m (see Ap- 
pendix[A]for the relevant NFW-equations), and then use Newton- 
Gauss Seidel relaxation with multigrid acceleration to obtain the 
solution. 

In spherical symmetry the field equation becomes one¬ 
dimensional, cj) = 0(r) and V 2 </> = ^4 + 2 For the Symmetron 
model, we implemented a simple fixed-grid multigrid solver for 
this one-dimensional problem. To do this we write Eq. as 
X/ = 0 where 


Xi = (VV); 


a 2 a ib Pm(n) 

2A 2 V a? pP 


\ 

- Vi + <p) 

/ 


(32) 


4 This condition is derived under the assumption of a spherical top-hat 
overdensity. 
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where subscript i denotes the value at the z’th gridpoint and 


(VV); = 


Vi +1 + Vi -1 “ 2 ipj 
(A r) 2 


2 Vi— 1 
r, A r 


(33) 


is a second-order discretization of the Laplacian. We start by 
making a guess for the solution, and then we loop through the 
grid updating the solution using 


Vi Vi ~ 


-Ci 

d-Ci/dipt' 


(34) 


The multigrid technique is used to speed up convergence. For 
more details about the methods we have used see iLlinares et al.l 
f2014} . 

The f(R) equation Eq. (13} is stiffer than the Symmetron 
equation for low values of /ro\, and a simple fixed-grid one¬ 
dimensional solver does not converge for the whole range of 
considered parameters. The main reason for this is that |/«| is 
required to be strictly larger than zero (but can get very close). 
This means that one has to ensure that |/r| does not become neg¬ 
ative in the solving process. To achieve this, we used the already 
well-tested /'(//(-module in the ISIS code ( |Llinares et al.j2014} . 
The only modification to the module is the change to a purely 
geometrical refinement criterion based on the value of p m . In to¬ 
tal, we used 10-15 levels of refinements starting from a base 
grid of N - 64 3 grid nodes in order to get accurate field profiles 
down to very small radii. The Symmetron is also implemented in 
the ISIS code, and we tested that our Symmetron field profiles 
agree very well with those obtained with our code. 

After having calculated the field profiles, the fifth force can 
be calculated by using Eqs. m and ( [12} for the Symmetron 
and f(R) model, respectively. For nDGP the analytical solutions 
to the field equation and force profile in spherical symmetry is 
given by Eqs. ([19} and ([20}, respectively. This makes it needless 
to solve any differential equations numerical I \f| The derivation 
of the analytical solution is shown in Appendix[B] 


3. Results 

In this section we present our results from numerical solutions 
of the field equations (§ |3.1| § |3.2[ ) and the obtained rescaling 
method for the three example models. The rescaling is first found 
empirically in §3.3 and then semi-analytically in the [mi a), [Hi a)) 
formulations in §[3.4| 


3. 1. Force profiles 

Figure [I] shows the radial profiles of y obtained solving 
Eqs. «[9f: 13 20 1 for the Symmetron, f(R), and, nDGP model, 
respectively (from left to right panels). We display a set of halo 
masses for each model using fixed model parameters. 

The differences in shape between the models are striking. 
The Symmetron model shows a whole range of solutions: from 
fully screened in the center for the heaviest clusters to not 
screened at all for the lightest ones. The maximum of the fifth 
force seems to be located ~ 1 — 10/Goo and drops again for larger 
radii. This drop happens due to the fifth force having a finite 
range, oc \e~' nr . 


5 Generally, this is true for all models that have a derivative shift sym¬ 
metry (like Galileons) and where the field equation is second order. If 
this is satisfied, then the equation of motion can be integrated to yield 
an algebraic equation in oc 


The f(R) model presented in the central panel of Fig. [T] in¬ 
herits a much more abrupt change from fully screened to not 
screened in the central halo region. With this particular choice of 
parameters (|/«ol = 1CT 6 * * , n = 1), halos with a mass of > 10 12 M Q 
seem to be fully screened in the central (R < 0.1 /Goo) region, 
whereas in lighter halos, the fifth force stays at its theoretical 
maximum 1/3 F^. For comparison, in Fig.[2]we show the force 
profiles derived using the semi-analytical method described in 
§ |3.4| for the f(R) (left) and Symmetron (right panel) models. 
The agreement is fairly good around the most significant region 
R ~ 1 - 5/Goo where y peaks. The semi-analytic results are pre¬ 
sented in more detail in § |3.4| 

In the righthand panel of Fig. [T] we show the radial scalar 
field profile tp(r) for a NFW halo in the nDGP model with r c = 1. 
As is also possible to see from Eq. <H]2}, the solution includes 
a curious feature: owing to the concentration-mass relation used 
(Eq. ( |A.3| )), heavier halos are less screened in the central region 
of the halo. Generally, the fifth force is less strong in the nDGP 
model and follows a similar functional form that is almost inde¬ 
pendent of mas^] 


3.2. Profile variation with halo mass 


To capture the change of the fifth force profile with halo mass, 
we chose to show the mass-weighted mean of the fifth force 


r 




dxy(x)w(x) 


(35) 


where w(x) is the normalized mass fraction w(x) d.r = d M/M(< 
x cut ) as defined in Eq. ( |A.4| >. This quantity was introduced by 
|Schmidt (|2010|l for the f(R) and DGP model and also considered 


Clampitt et al. (2012} for the symmetron model. The outer 
cutoff was fixed to x cut = 10. This value was chosen so that 
vv(Xcut) is roughly one-tenth of the maximum weight, meaning 
that not too much information is lost. As discussed in § |4.2| (see 
Eq. ([55}), y can bea related to the ratio of the dynamical mass to 
the lensing mass. 

Figure [3] shows y versus the halo mass. However, this mea¬ 
sure is not unique. Instead it is also possible, for example, to use 
a difference between the force in the center and in the outskirts 
of the halo. We discuss the advantages and caveats of using y as 
an “observable” in § |4.2| 

The lefthand panel of Fig. [3] shows three Symmetron solu¬ 
tions with different model parameters (a ss b, L, fi). The slope, as 
well as the offset, of the three curves is seen to be fairly dif¬ 
ferent. Also notable is the immense change in y with changing 
halo mass. The f(R) solutions (shown in the central panel of 
Fig. 0 possess an upper limit at the already mentioned value of 
one-third. For lighter halos y falls off rapidly. In the righter-most 
panel of Fig. [3] we show the variation in y with mass for the 
nDGP model for three different values of r c . Overall, the impact 
of the fifth force is much weaker than in the other models. As 
discussed in m the variation in the fifth force with halo mass 
is only due to the mass-concentration relation of the NFW halo. 
This feature manifests itself in Fig.[3]in the slight increase in y 
with increasing halo mass. 


3.3. Empirical rescaling 

To find a rescaling that maps the various characteristics pre¬ 
sented in the last section onto one single curve, we need to 

6 The only mass dependence is implicit in that NFW halos of different 

sizes have different values for the concentration c(M). 
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Fig. 1. y vs R/R 200 for NFW halos with masses ranging from 10 9 M a h 1 to 10 14 ' 5 M a h 1 . The left panel shows the Symmetron solution using 
(a ss b, L, /3) = (0.7, 1 Mpc/z -1 , 1), the central panel the Hu & Sawicki f(R) solution with (|/ fl0 ], n ) = (10 6 , 1), and the right panel the nDGP 
solution with r r = 1. 


Symmetron 


/(«) 



R/ R200 


Fig. 2. Semi-analytical solutions for y for the Symmetron (left) and f(R) 

rescale the x and the y axes of Fig. [3] The latter is naturally nor¬ 
malized to the theoretical possible maximum of the gravitational 
enhancement y max . The rescaling of the halo mass is less triv¬ 
ial. Therefore, we simply define a characteristic mass p^m in the 
following way: 

7<>200) = ^7max (36) 

where y is still the mass-weighted average of the enhancement of 
the fifth force as defined in Eq. ( [35] ) In general ^200 is a function 
of the various model parameters. Halos with mass much higher 
than / 22 oo are in the fully screened regime, whereas halos with 
mass much lower than 1/200 are in the unscreened regime. 

In this section we find the functional form of yi /200 for each 
model merely empirically from our numerical solutions. How¬ 
ever, in the next section, we derive semi-analytic predictions for 
a large class of scalar tensor theories. 

This means that, by looking at Fig. [3] it is possible to read 
off several values of // 200 - For example, for the f(R) model with 
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R/R 200 


(right) using the same parameter values (and color coding) as in Fig. [T] 

(n, l/sol) = (1, 1(T 6 ) (in green), y(~ 10 13 M 0 /r‘) = 1/6 = 
1 / 2 y max , i.e., 7 / 200 . is in that particular case ~ 10 u M Q h 1 . 

For the Symmetron model, we found //200 to be well fit by 

A” = 2 X 10 , 0 Mo*-' X J . 07) 

and y max is given by Eq. GD- Figure [4] shows the resulting 
rescaled curves, i.e., y/y m!a versus M/// 200 • In addition to the 
four parameter sets already presented in the lefthand panel of 
Fig. [3] (keeping the color coding), we plot an additional 125 
curves with randomly chosen parameters. Specifically, we drew 
L ~ [0.2, 2]Mpc hr 1 , a ssb ~ [0.2, 0.8], and p ~ [0.2, 2.4] (all 
uniformly). 

It is notable that in spite of this wide parameter range, the 
resulting curves in Fig. [4] are very similar. All follow the char¬ 
acteristic shape - from y ~ y ma x at M <sc //200 through y(M ~ 
// 200 ) = 7max/2 to y ~ 0 for M » P 200 - with very little scatter. 
We note, however, that for low masses, y does not fully approach 
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Symmetron f(R) nDGP 



Fig. 3. y vs halo mass M. The left panel shows the Symmetron curves with (a ss b, L , (3) being (from solid to more interceptions) 0.7, 1 Mpc h l , 1 
(red),0.7, 2Mpc/j _1 , 1 (blue), 0.3, lMpc/;' 1 , 1 (green), 0.7, IMpc/? -1 , 2 (purple). The central panel shows the;Hu & Sawicki f(R) solution with 
log , 0 l/fiol = (-4, -5, - 6 , -7) (in red, blue, green, and purple, respectively) and the right panel the nDGP solution with r c = (0.5, 1,3) (from top 
to bottom). 


y max but falls ~ 10% short. This is due to the drop in y for larger 
radii as is visible in the lefthand panel of Fig. |T] 

In Fig. [5] we present the results of the same rescaling tech¬ 
nique applied to the f(R) results. In this case the rescaling mass 
is found to be 

^ ) = 10 13 M o r 1 x(Mj (38) 

and y max = 1/3. The curves in Fig. [5] correspond to the numer¬ 
ical solutions for each combination of n = ( 1 , 2 ) and |/ R0 | = 
(10~ 8 , 3 x 10~ 8 , 10- 7 , 3 x 10- 7 ,..., 1(T 4 ), with the ones already 
presented in Fig.[3]color-coded accordingly. Also here, the scat¬ 
ter around the mean form is quite small. The second striking 
feature the similarity in the overall shape of the curves in Figs. [4] 
and [5] — even though the solution stems from two independent 
modified gravity theories. The main differences between the two 
models are the slightly different transitions widths and, the fact 
that for M <s // 200 , 7 ~ Tmax in the f(R ) case, whereas the Sym¬ 
metron curves are limited to a lower value. 

Rescaling the 7 -mass relation for the nDGP model (pre¬ 
sented in the right panel of Fig. [3} is not as straightforward. Since 
the curves are nearly constant (and do not possess a characteristic 
drop), no empirical formulation for 7/200 is apparent. However, 
since calculating 7 in the nDGP model is not computationally 
expensive, we can numerically solve Eq. ( |36| ) for several values 
of r c . This let us fit the empirically found relationship 

logioTDoo = A - (r c /a)V' v/a , (39) 

which yields (A, a, b) = (11.14, 6.72, -0.57). The rescaled 
y(M 2 oo) solutions obtained with r c = (0.5, 1, ..., 20) showed so 
little variation that we chose not to explicitly show them. Instead, 
we display in Fig.[ 6 ]the average y/jmnx - M 200 /R 200 curves (solid 
lines), as well as the maximal deviation found for each (dashed 
lines). Clearly, the nDGP curves show very little variation. 

Figure [ 6 ] once more shows the difference between the Sym¬ 
metron and f(R) model with their characteristic ‘inverse-S- 
shape’ on the one side, and, the nDGP - with a log-linear relation 



Fig. 4. 7 versus M/p 2 m f° r the Symmetron model. The solid colored 
lines indicate the same model parameters as in Fig. [3] The gray lines 
are an additional 124 randomly chosen Symmetron parameters. 7/200 f° r 
the Symmetron is given by Eq. ( |37| ). The black dashed lines indicate the 
point defining p 2 oo in Eq. ( |36| ), and the gray arrows point to the screened 
and unscreened regimes. 


- on the other. A simple fitting function that is able to match the 
characteristic ‘inverse-S-shape’ given by 


y(M) = y max 


1-arctan(M 2 oo/// 20 o)]. 

71 


(40) 


However, it is also clear that the models of the former group 
do not show exactly identical behavior. Instead, the transition 
region for the Symmetron model is a bit wider, and, as mentioned 
before, it does not converge to y max for lower masses. This is why 
we propose the more general form 

y(M) = y max (1 - 


- arctan [(M 2 oo/R 20 o)° 

71 


(41) 


with two free parameters JR, a. Here, the prefactor JR can be in¬ 
terpreted as the part of the theoretical enhancement of gravity 
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Fig. 5. y versus M/p 2 oo for the Hu-Sawicki f(R) model. The solid col¬ 
ored lines indicate the same model parameters as in Fig. [3] The gray 
lines are an additional 5 values of f R 0 . ^200 for this model is given by 

Eq. {38}. 



Fig. 7. Semi-analytical predictions for y versus M/p 200 f° r the Hu- 
Sawicki f(R ) model (red), the Symmetron (blue) and the Dilaton model 
(purple). The solid black curve shows the fitting formula Eq. {40}. The 
shaded regions show the results of the full numerical analysis. 



Fig. 6. y versus M/p 200 f° r the three example models. The solid 
[dashed] lines show the mean [minimum & maximum] values for the 
128 (Symmetron), 16 (/((?)), and 20 (nDGP) numerical curves ob¬ 
tained. The solid gray lines show the fits of Eq. <ED described in jj3.3| 


that is effectively active (and, thus, can only be < 1), and a con¬ 
trols the width of the partially screened region with a = 0 being 
‘infinitely wide’ and a —> 00 the transition turning into a step 
function. Also, if a < 0, the “inverse-S-shape” flips vertically 
and turns into an “S-shape”. Quantitatively, we can relate a to 
the half width of the partially screened region through 


W(a) = tan j . (42) 

Here, M 200 € Ipzoo/W, Wp 200 ] is the interval where d < 
r/Tmax <3\(l -d). 

Figure [6] also shows the results of fitting Eq. {41} to all of 
our numerical data (with enforcing JK < 1). It is noteworthy that 
this simple, two-parameter fit describes the numerical results of 
the three studied models very well, even for the nDGP model. 
As expected, the a parameter reflects the fundamental difference 
between the nDGP and the two other models, as it is close to zero 
(-0.03) for the nDGP and more similar for the other two models. 
Specifically, we obtained a = 0.46 and 0.78 for the Symmetron 
and f(R) model, respectively. 
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3.4. Semi-analytical derivation of p 200 in the { m(a),j3(a )} 
formulation 


In the following, we discuss the origin of the empirical rescalings 
found. Using the general framework of the {J3(a ), m(a)\ formula¬ 
tion described in j {2.3| we can derive an analytic approximation 
for p 2 oo- It turns out that we can do this quite generally, so our 
results below can be applied to any scalar-tensor theory with a 
Lagrangian L = \{dyy L + V(ip) and a conformal coupling to mat¬ 
ter. We only focus on a very rough derivation, which can be seen 
as the foundation of a more precise treatment in the future. 

Before we do this, we should check to what extent our semi- 
analytic approximation, which formally holds only for a spher¬ 
ical top-hat overdensity, is valid for NFW halos. In Fig. [2] we 
show the semi-analytical force profiles that can be compared to 
Fig. 0 which shows the true numerical results. The qualitative 
agreement is fairly good. We are able to match the shape around 
the peak, which is the significant part of the profile that will give 
rise to most of the signal in y, and the amplitude is off by no more 
than a factor of a few, which is good enough for our purposes. 

As a first - rather crude - approximation, we assume that a 
halo is fully screened within a certain radius x scr een. In the un¬ 
screened regime x scr een —> 0, but generally 0 < x scr een $ 10. 
Thus, we approximate y(x) ^ ymax#(*screen - x) where 8 is 
the Heaviside function. This leads to a solution for the mass- 
weighted average of y, which reads as 


y ~ 


Tmax 


A7( < X s creen) 

47(< x cut ) 


(43) 


Furthermore, using the definition of P 20 Q, namely y{p 2 ao) = 
57 max, we require 


M (< -^screen) ^ 1 

M{< X cu t) 

1 ^ 200—^200 ^ 


(44) 


which requires x scr een ~ (9(1) for typical halo masses. 

The dependence of x scrC en, hence of p 2 oo, on the various SMG 
model parameters can be calculated using the screening condi¬ 
tion Eq. {30}. In this case, the condition reads as 


\(fi(a p (x)) - ip(a 

env)l 

2/?(#env )^P1 (-^screen) 


^p(^) — ^env 


- \ 1/3 

Pm 


p(x 

screen) 


(45) 

(46) 
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With Oiv(r) = GM ^ <I 2 and M 200 = 47r/3A v i r Q m p c or 3 —> r oc M 1 ^ 3 , 
we can easily calculate how /i 200 scales with the model parame¬ 
ters in any given model. 

For the Symmetron we find the scaling M(< x scr een) 2 ^ 3 ^ t - - 

a ssb 

As a result, P 200 scales as 


Symmetron _ , 

^200 _ ' 


0.4 ± 0.3) x 10 10 Ma/r 1 x 


ya^Mpc/h 


(47) 


where the proportionality factor was calculated numerically. 
This scaling is exactly the same as we found empirically in the 
previous section, but the proportionality factor is off by a factor 
of-3-10 . 

For the Hu-Sawicki f(R) model, we find the scaling P 200 
l/flol 3 ^ 2 C 3/,2< " +1) fW 0 //!, where C is a constant of order one, depen¬ 
dent only on x screen and G ,„■ Numerically we find that the results 
are only weakly dependent on tv. changing n from 1 to 5 only 
gives us a factor 1.5 offset in the value we find for // 200 - Again 
the proportionality factor is found numerically and gives us 


/4oo' = (~ 11 ± 0-4) x 10 13 M Q r' 


x 


mr 


(48) 


which is almost spot on the empirical relation found in the pre¬ 
vious section. 

For a general {m(a),/3(a)} model, we see from Eq. ( [23] ) that 
we would (again very roughly) expect to have (f>{a) 


m 


a 3 m 2 (a) * 

This in the condition Eq. © means that we generally expect 
P 200 K 4^(a) where A<p{a) = rrT l {a) is the range of the fifth force. 

In Fig.[7]we show the semi-analytical profiles for y in terms 
of M 200 //U 00 for f(R) and the symmetron, together with the full 
numerical results. To demonstrate that this rescaling works in 
general, and not just for the two models presented here, we also 
calculated the y predictions for the Dilaton model presented in 
Brax et al. (2012c |2010) >. This model is characterized by 

m(a ) = moa~ r (49) 

P(a)=/3 0 e^ 2r - 3 ^) (50) 


where r, s, mo, and f3o are free model parameters and where the 
evolution of /3 are very different from our two example models. 
Figure [7] shows that the rescaling also works perfectly for the 
Dilaton model. 

We have not been able to get a prediction for the width of the 
partially screened region semi-analytically, so we leave this for 
future work. 


3.5. Predictions for other screening mechanisms 


For other screening mechanisms that are not encapsulated by 
the { tn(a),/3(a )) formulation, it is harder to make general pre¬ 
dictions. For any theory with a derivative shift symmetry (the 
Galileon symmetry), second-order equations of motion the field 
equation in spherical symmetry are integrable, leading to an al¬ 
gebraic equation on the form 





A vir 

/logd+CX)-^ 

3x 3 

l logd+c)-^ 


(51) 


where the last fraction can be identified as gNFw(cx)/gNFw(x) 
(see Eq. (|A.2|>). For a general model this leads to 


7 = Tma xG 


gNFw(c) 

X 3 gNFw(cx) 


(52) 



9 . 10 . 11 . 12 . 13 . 14 . 15 . 


Log 10 (M h/M o) 

Fig. 8. y(M) for NFW halos with masses ranging from 10 9 M q /T 1 to 
10 15 M a tr l using Eq. {52} with G{z) = 2((1+ ^ > — . The four bands of 
lines show the results (from top to bottom) of £ = 0.1,1.0,10,100, and 
within each band, we have (from top to bottom) n = 0.5,0.4,0.3,and 
0 . 2 . 


for some function G that satisfies G(0) = 1 and GY 00 ) = 0. For 


DGP we have G(z) = 


2 (Vl+fz-D 


with £ = 


8 cH;n,„A,j 


Y ^ 9/^DGP 

For models that only have a shift symmetry, like the k- 
mouflage screening mechanism, the field equation. 


1 d 


--[r z /( 2 Of] 


d(f) PSp n 


dr 


dr 


Mp\ 


(53) 


can also be integrated up, and from this it follows that the fifth- 

r 8 m r 2 dr 

force /-j 4 oc V0 is a function of ;iL -p— for spherical symmetry. 
This leads to a slightly different form 


T = Tma xF 


I gNFw(c) ) 

\X 2 gNFw(cx)/ 


(54) 


for some function F satisfying F(oo) = 0. 

The exact form for the functions F and G are model depen¬ 
dent, making it hard to make definite statements here. However, 
one general feature we are able to deduce is that, just as we saw 
for DGP, y is only sensitive to the concentration of a halo. We 
therefore expect a fairly weak mass dependence of y for models 
in this class. 

As an example, we can try to parametrize the functional form 
for F and G and calculate y. One functional form that satisfies 
the requirements is G(z ) = (1+ (; 3 1 with n e (0,1). The results 
we get from this of course depend on the values we take for 
the parameters (like £); however, taking realistic values for the 
parameters, which is motivated by the requirement of satisfying 
local gravity constraints, the typical result we find is as expected 
(see Fig. [8]>: a very weak M dependence of y. 


4. Discussion 


Our numerical force profiles presented in Sects. 0 and |3.2 


are 

in good agreement to what has been previously found in A'-body 
simulations (e.g., |Zhao et al.|201 la[|Falck et al.|2015} . The sim¬ 
ple rescaling procedure based on the maximum enhancement of 
the gravitational force y max and the intermediately screened mass 
AUoo ( §3.3} demonstrates that in spite of their foundational differ¬ 
ences, the three screened-modified gravity models investigated 
show self-similar behavior in cluster environments. 
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Other parameters: 



Fig. 9. Rescaling parameter P 200 versus the maximum gravitational en¬ 
hancement y max for several discussed models (solid colored lines) and 
the respective partially screened regions (colored semi-transparent ar¬ 
eas). The red line shows the space spanned by the Hu-Sawicki f(R) 
model when varying l^o] from 10~ 8 (left end) to 1(T 4 (right end), while 
leaving n = I fixed. For the Symmetron, we also varied one parameter 
per curve (in blue, green, and violet), while leaving the others fixed at 
= 1, L = IMpcfi -1 , a ssb = 1/2. The orange line indicates the region 
for nDGP models with r c e [0.5, 10] (left to right) with the associated 
semi-transparent area spanning wider than the plot range. See j ]4.1| for 
details. 

4.1. Observational implications 

Our rescaling methods allow an abstract characterization of 
screened modified gravity theories, which themselves might pos¬ 
sess quite different screening mechanisms. This potentially has 
interesting implications for constraining SMGs. 

The colored lines in Fig. [9] show y max versus P 200 for some 
selected models. The lines were drawn by altering one model 
parameter while leaving the others constant (see figure cap¬ 
tion for details). Each set of SMG model parameters occu¬ 
pies a specific point on this graph; for example, the combina¬ 
tion (a ssb , P, L ) = (1/2, 1, IMpc/r 1 ) maps to (y max , P 200 ) 
(1.75, 4.5 x 10 n M G /! _1 ). The semi-transparent areas in Fi g. |h| 
denote the extent of the partially screened region using Eq. \42) 
with d = 0.2. This means that for every selected value of /./ 200 , the 
semi-transparent area indicates the range where the fifth force is 
at least 20% but a maximum of 80% active. 

Thus, Fig. [9] indicates (i) what halos are partially screened 
and, (ii) what the maximum enhancement of gravity is for this 
model. Physically, P 200 /W can be interpreted as the minimum 
mass where one expects a variation in halo properties; that is to 
say: 

- For M 200 £ P 200 /W halos in this model are unscreened, 
thus allowing for comparison with the ACDM predictions. 
In particular, one expects that the majority of halo proper¬ 
ties are rescaled by ~ y max . This can make the comparison 
with ACDM difficult when using, for example, the dynami¬ 
cal mass as mass estimate. 

- Halos with M 200 e [P 200 /W, Wp 200 ] are partially screened. 
This means that in this mass range, halo properties might 
differ from the expected. Since halos in this mass range are 
between screened and unscreened, the environment of these 
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halos has a major effect (see|Zhao et al.|201 lb 1 Winther et al. 


2012 for the study of this effect using /V-body simulations). 


Also, comparing the inner and outer regions of individual ob¬ 
jects can lead to insight (for instance, the ratio Msoo/A^oo). 
In Fig. [9] the partially screened region is shown as a shaded 
area in the matching color. As demonstrated for the nDGP 
case, it is possible that - for some SMGs - all halos are al¬ 
ways partially screened. 

- For masses > Wp 200 - the halo is fully screened, so no devia¬ 
tion from the ACDM predictions is expected. 

Conclusively, in any observable-mass scaling relation, we pre¬ 
dict a peculiar feature at P 200 , where the observables match the 
ACDM prediction for higher masses and where they are differ¬ 
ent for lower masses. However, the latter statement only holds 
if the true mass of the halo is used. This is because mass tracers 
like the dynamical mass can also be affected by the fifth force 
dSchmidt|2010| >. 

The extent of the partially screened region, which we quanti¬ 
fied with a (or, W , see Eq. ( |42| , can potentially be used to distin¬ 
guish between various screening mechanisms. We found that this 
extent is constant for each screened-modified gravity considered, 
but differs significantly between them. In particular, the value of 
a found in the numerical fits of 03] translates to log 10 W ~ 16, 
~ 1, and ~ 0.6 for the nDGP, Symmetron, and the f(R) models, 
respectively (using d = 0.2, as above). In a y max -|oj (or equiva¬ 
lently, a y max -W) diagram, we expect eventual constraints to rule 

OUt the area y max > yconstraint' kr 5) ^constraint (W > W(ctconstraint))' 


4.2. Caveats 

Several assumptions went into deriving the presented formalism. 

Here, we want to mention and investigate them. 

- Halo profiles and shape. Throughout this work, we assumed 
spherical symmetry and perfect NFW halos. However, in re¬ 
ality clusters of galaxies do not exactly follow an NFW mass 
distribution, and halos also possess a nonzero ellipticity (e.g., 
|Oguri et al.|2010| >. Another phenomenon falling in the same 
category is our assignment of a unique mass-concentration 
mapping, whereas the spread in this relation is fairly wide 
(Neto et al. 2007) i. Both are consequences of the fact that 
clusters are individual objects with their own histories and 
environments. As a result, individual clusters (as any astro- 
physical object) have mostly a limited explanatory power and 
must be stacked for analysis. 

- Model assumptions. To rescale the y - /1200 curves accord¬ 
ingly, an intermediate-mass scale has to exist where y ~ 
ymax/2- Since this is the scale between the screened and un¬ 
screened regimes and, our work addresses ‘screened mod¬ 
ified gravity theories’, we think that this can be safely as¬ 
sumed. The universal rescaling works perfectly for the mod¬ 
els we have tested, and we conjecture that it will hold for 
all { 111 (a),/3(a)] models. However, this might not be true for 
cases where m(a ) and/or /3(a) are very steep functions of a, 
leading to a very rapid evolution of the screening condition 
with density. 

- Choice of parameters. Throughout this work, we fixed a 
number of parameters. Probably, the most influential ones 
are the halo cutoff radius x cut = 10 and defined P 200 as 
y(P 20 o)/yrnax = k with k = 1/2. The latter choice is arguably 
quite natural. In spite of that, we do not expect our results to 
change dramatically if a (slightly) higher or lower value of k 
is taken since the model curves are expected to shift equally. 
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The choice of x cut affects our results more drastically, as can 
be inferred from Fig. |T] A lower value of x cut leads, for in¬ 
stance, to a decrease in the extent of the partially screened 
region (i.e., an increase in a). Nevertheless, since all models 
are affected by this, the precise choice of x cut (within reason) 
does not alter our conclusions. 

- Practicality. We focused merely on the force profiles di¬ 
rectly, leaving the question how practical the study is since 
y(r) and y are not directly observable. Although this is not 
wrong, we caution that the fifth force triggers many observ¬ 
able phenomena, as has been shown in various /V-body sim¬ 
ulations. In particular, y(r) can be related to the ratio of the 
dynamical mass to the lensing mass (|Schmidt||2010f |Zhao| 
|et al.|201 lbl|Winther et al.|2012| > 


y(j) = 


M D yn(< r) 
^Lensing (< Q 


(55) 


where 



dr- 

(56) 

dr. 

(57) 


Here, (® + 'Fj/2 is the dynamical potential, (® - M J )/2 is the 
lensing potential, and ® and 'F are the two metric potentials. 
For the models considered in this paper, the lensing poten¬ 
tial is the same as in ACDM, but there are several models 
where the lensing potential is modified as in the Galileon, 
for example^] In practice, however, it is not possible to mea¬ 
sure the potentials directly. Instead, the dynamical mass is 
usually obtained via the cluster velocity dispersion through 
a power-law relation obtained from (ACDM) /V-body simu¬ 
lations (e.g., |Becker et al.||2007| ). The lensing mass can be 
read from the stacked density profiles assuming a constant 
mass-richness relation ( Johnston et al.|2007||Hoekstra et al.| 
|2013| >. This means that our mass-weighted average of the 
gravitational enhancement y is directly proportional to the 
ratio of the dynamical to the lensing mass as has been shown 
by 


Schmidt 


(2010 


In conclusion, our work can be seen as a first step toward a truly 
unified description of SMGs on cluster scales. Naturally, this first 
step can be tuned further and extended to improve its precision 
and include more models. 


5. Conclusions 

Starting from the simple idea that screened modified gravity the¬ 
ories in general possess three regimes within nonlinear astro- 
physical structures, a fully screened, an unscreened, and a par¬ 
tially screened regime, we investigated whether one can spec¬ 
ify, independently of the SMG model, at which halo mass range 
these regimes are located and how extended they are. 

These two questions can be answered as follows: 


- It is possible to formulate the mass of the partially screened 
regime as a function of various model parameters. We 


7 While this paper was in the final stages of completion, Barreira et al. 
(2015 ft showed that for the Cubic Galileon model, the lensing potential 
within clusters is practically unmodified. 


In fact, Schmidt 


l) 3/5 M L en S mg (see their Eq. (64).). 


(2010' derive the relation M Dyn ,observed = (y + 


demonstrated this empirically for the Symmetron, the Hu- 
Sawicki f(R), and the nDGP models ( §3.3| >. In addition, we 
found a semi-analytic derivation of this functional form in 
the \m(a), /J(a))-parametrization of Brax et al. (2012b I. This 
allowed us to expand our method easily to the Dilaton model. 
The extent of the partially screened regime is an intrinsic 
property of a particular screening mechanism. We found that 
the mass range of the partially screened region differs signif¬ 
icantly between the three example models but are (approxi¬ 
mately) constant within a particular model. 


In spite of the differences in the nature of their screening 
mechanisms (see § |2. 1 1 ), we characterize the Symmetron, the Hu- 
Sawicki f(R), and the nDGP models with three parameters: (i) 
the center of the partially screened regime P 200 > («) the maximal 
enhancement of the fifth force y max , and the (Hi) effectiveness 
of the screening a (which controls the extent of the partially 
screened regime). We derived the first two parameters analyt¬ 
ically in the {/«(«), /ii«))-framework, which includes the Sym¬ 
metron, the Hu-Sawicki f(R), and the Dilaton models. 

Quantifying observational constraints in the proposed 
(P 200 > Tmax, aO-parametrization thus allows several screened 
modified gravity models to be constrained in bulk instead of hav¬ 
ing to consider them - and their respective parameter space - 
individually. 
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Appendix A: NFW equations 


The Navarro-Frenk-White (NFW) profile is an empirically de¬ 
rived halo profile fNavarro et al.|1995 [Navarro etal.|1996l|1997| l 


from /V-body simulations. The NFW density profile reads as a 
function of normalized distance from the center x = r/R va : 

p(x) _ Ayi r c 2 g NF w(c) 

Pm 3x( 1 + ex) 2 
where 


(A.l) 


£nfw(c) - 


1 


ln(l +c)-c/(l +C) 


(A.2) 


Here, c is the concentration parameter, which we link to M 200 
through the observed relationship found in|Bullock et al.|([2001|>: 


c = 9 


M200 


3.2 x 10 l2 M a h- 1 


- 0.13 


(A.3) 


Throughout this work we define A vlr = 200, which means R vll = 
R 200 and M vir = M 20 o = 200 x 4/3 xR^PcO- 

Although the density of the NFW profile is diverging for x —> 
0, the enclosed mass within a certain radius does not (jLokas &| 
|M amon |200T| ). This leaves the mass fraction within [x ± dx/2] 
as 

d M _ c 2 xgNFw(cx cut ) 

M« x cul ) “ (1 +cx) 2 ’ ( ' J 

where x cut is defined below Eq. ( |35j ). This will be useful when 
considering mass-weighted quantities. 

To calculate the screening condition we also need the expres¬ 
sion for <I> v. For a NFW halo, we find 

IrT. 1 ( \ (^0^vir)~A v j r gNpw((') 1 /1 , /. c\ 

l®jvl(x) =-log(l + cx) . (A.5) 
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Appendix B: DGP solutions for NFW halos 

In this appendix we present the equations for y for the DGP 
model. 

For a NFW density profile (Eq. ( |A. 1 1 >) halo, we have 

Id 5 <n(r)r 2 dr _ Avu - / iog(l + cx) - _ A vir g NFW (c) 

r 3 3x 3 \ log( 1 + c) - / ^ 3 g NF w(cx) 

where x = j-. This leads to the fifth force profile being given 
by 


7 = 7ir 


2 ( Vl + e(x) - 1) 


e(x) 

where y max = ^ and 

8 (r f // 0 ) 2 G„, A vir g NFW (c) 


e(x) = 


9/^dgp * 3 £nfw(cx) ' 


(B.2) 


(B.3) 
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